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\ Abstract 

We consider the problem of sparse estimation in a factor analysis model. A traditional 
estimation procedure in use is the following two-step approach: the model is estimated 
by maximum likelihood method and then a rotation technique is utilized to find sparse 
factor loadings. However, the maximum likelihood estimates cannot be obtained when the 
number of variables is much larger than the number of observations. Furthermore, even if the 
maximum likelihood estimates are available, the rotation technique does not often produce 
a sufficiently sparse solution. In order to handle these problems, this paper introduces a 
penalized likelihood procedure that imposes a nonconvex penalty on the factor loadings. 
We show that the penalized likelihood procedure can be viewed as a generalization of the 
traditional two-step approach, and the proposed methodology can produce sparser solutions 
' than the rotation technique. A new algorithm via the EM algorithm along with coordinate 

o ' 

£NJ ■ descent is introduced to compute the entire solution path, which permits the application to 

a wide variety of convex and nonconvex penalties. Monte Carlo simulations are conducted 
to investigate the performance of our modeling strategy. A real data example is also given 
C3 , to illustrate our procedure. 

Keywords: Coordinate Descent Algorithm, Nonconvex Penalty, Rotation Technique 



1 Introduction 

Factor analysis provides a useful tool for exploring the covariance structure among a set of ob- 
servable variables by construction of a smaller number of unobservable variables. In exploratory 
factor analysis, a traditional estimation procedure in use is the following two-step approach: the 
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model is estimated by the maximum likelihood method wit h the use of efficient algorithms (e.g. 



Joreskogl . 



1967 



Jennrich and Ro 



as the varimax method (IKaiser 



bmson, 



1969 



Clarke 



19701 ). an d then rotation techniques, suc h 



19581 ) and the promax method (IHendrickson and Whitd . Il964j ). 
are utilized to find sparse factor loadings. However, it is well known that the m aximum likeli hood 
method often yields unstable estimates because of overparametrization (e.g., 



Akaike, 



19871). In 



particular, the above-mentioned algorithms cannot often be used when the number of variables 
is much larger than the number of observations. Furthermore, even if the maximum likelihood 
estimates are available, the rotation technique does not often produce a sufficiently sparse solu- 
tion. In such pena lized likelihood procedure that produces the sparse solutions, such as 



the lasso (ITibs 



hiram 



The lasso ( ITibshiranil . 



1996 ]), can handle these issues. 



19961 ) has received much attention as a powerful tool for sparse es- 



timation in a variety of statistical modeling and machine learning, includin g generalizec 



models, Gaus s ian graphical models, and t he support vector machine (e.g. 



Hastie et al. 



mear 



2004 



Yuan and Linl . 



2007 



Friedman et al. 



2010 



and estimates an overly dense model ( 



Zoul . 



. It is , however, we 



2006 



Zhao and Yu 



1 known that the lasso is biased 



2007 



ZhaneJ . 120101 ). Typically, 



a penalization methods via nonconvex penalties can achieve sparser models than the lasso. Re- 
cently, a number of researc hers have presented various nonconvex methods: brid ge regression 



2001 



( Frank and Friedman! . 



1993 



Fu 



19981 ). s moothly clippe d absolute deviation (SCAD, 



minimax concave penalty (MC+, 



Zhang , 



Fan and Li . 



20101 ) . and generalized elastic net ( jFriedmanl . 



20121 ). Several efficien t algorithms to obt ain the entire solutions hav e also been propo sed (e.g., local 



linear approximation, IZou and Li 



ear unbiased selection a 



2011 



Mazumder et al 



gorit hm, 



2008; generalized path seeking, 



Friedman 



Zhang , 



2010 



coordinate descent algorithm, 



2012 



; penalized lin 



Breheny and Huangi . 



201 lh . 



In the framework of factor analysis m o dels, a few researchers 
likelihood procedure 



Ning and Georgioul ( 120111 ) and 



Choi et al 



rave discussed the penalized 



(j201l[ ) applied the lasso-type 



penalization procedure to obtain sparse factor loadings, and numerically showed that the penal- 
ization method often outperformed the above-mentioned two-step traditional technique. However, 
the relationship between the penalized likelihood procedure and the two-step traditional approach 
has not yet been discussed. Moreover, although the ordinary lasso can estimate an overly dense 
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model in penalized likelihood factor analysis (iChoi et al. 



2011 



Hirose and Konishi 



20121 ) as dis- 



cussed earlier, the penalized likelihood methods via the nonconvex penalties that produce sparser 
solutions than the lasso have not yet been developed. 

In the present paper, a new penalization method via nonconvex penalties in factor analysis 
models is proposed. We show that penalized likelihood method can be viewed as a general- 
ization of the traditional two-step approach, i.e., the maximum likelihood estimates with the 
rotation technique can be obtained by the penalized likelihood procedure under certain condi- 
tions. The proposed methodology can produce sparser solutions than the rotation technique 
with maximum likelihood estimates by changing the regularization parameter. In order to pro- 
duce entire sol utions of factor loadings a nd unique variances, a new pathwise algorithm via the 



EM algorithm ( 



Rubi n and Thayer 



( IMazumder et al 



1982j) along with coordinate descent for nonconvex penalties 



20111 ) is introduced, which permits the application to a wide variety of con- 
yex and nonconvex penalties inclu ding the lasso, SCAD, and MC+ family. A package fane in R 
(|R Development Core Teaml . 120081 ). which implements our algorithm, is available from Compre- 
hensive R Archive Network (CRAN) at htt p: //cran| . r-pro j ect . org/web/packages/f anc/index . html. 

The remainder of this paper is organized as follows: Section 2 describes the traditional estima- 
tion procedure in factor analysis. In Section 3, we introduce a penalized likelihood factor analysis 
via nonconvex penalties, and show that penalized likelihood procedure can be viewed as a general- 
ization of the maximum likelihood method with the rotation technique. Section 4 provides a new 
algorithm based on the EM algorithm and coordinate descent to obtain the entire solution path. 
In Section 5, we present numerical results for both artificial and real datasets. Some concluding 
remarks are given in Section 6. 



2 Traditional Estimation Procedure in Factor Analysis 



2.1 Maximum Likelihood Factor Analysis 



Let X = (Jfi, . . . ,X p ) be a p-dimensional observable rando m vector with mean vector fx and 



variance-covariance matrix S. The factor analysis model (e.g., 



Mulaik 



20101 ) is 



X = fi + AF + e, 



where A = (Ay) isapxm matrix of factor loadings, and F = (Fi, ■ ■ ■ , F m ) T and e = (ei, ■ ■ ■ , e p ) T 
are unobservable random vectors. The elements of F and e are called common factors and unique 
factors, respectively. It is assumed that the common factors F and the unique factors e are 
multivariate-normally distributed with E(F) = 0, E(s) = 0, E(FF T ) = I m , E(es T ) = and 
are independent (i.e., E(Fe T ) = O), where I m is the identity matrix of order m, and VP is a 
p x p diagonal matrix with z-th diagonal element ipi, which is called unique variance. Under 
these assumptions, the observable random vector X is multivariate-normally distributed with 
variance-covariance matrix X = AA T + ty. 

Suppose that we have a random sample of N observations X\, ■ • ■ ,xn from the p-dimensional 
normal population N p (fi, AA T + V&). The log-likelihood function is given by 



(A, *) = -— [plog(27r) + log | A A' 



T 



*| +tr{(AA T + *)- 1 5}] 



(1) 



where S = (s^) is the sample variance-covariance matrix. 



The maximum likelihood estimates of A and ^ are given as the solutions of d£(A, ty)/dA = 
and d£(A, ^) / d^f = 0. Since the solutions cannot be expressed in a closed form, several researchers 
have proposed iterative algorithms to obtain the maximum like l ihood estimates Aml and ^ml 



(e.g., 



Joreskog . 



19671 : 



Jennrich and Robinson! . 



1969 



ClarkeJ, 



197G; 



Rubin and Thayer 



19821 ). 



2.2 Rotation Techniques 

The factor loadings have a rotational indeterminacy, because both A and AT generate the same 
covariance matri x S, where T is an arbitrary orthogonal matrix. The rotation technique (e.g., 



varimax method. 



Kaiser 



19581 ) has been widely used to find the matrix T which gives a meaningful 



relation between items and factors. Suppose that Q(A) is an orthogonal rotation criterion at A. 
The criterion is minimized over all orthogonal rotations with an initial loading matrix being Aml, 
i.e., 



minQ(A), subject to A = A ML T and T T T 

A 



(2) 



For ease of comprehension, here and throughout this paper, the criterion Q (A) is given by 



the component loss criterion Q(A) = Y^i=i Sj=i -f(l^iil) (jJennrich 



2004 



20061 ). Note that it is 



not necessary to assume that Q(A) is a component loss criterion for constructing our modeling 
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procedure. The problem in ([2]) is then expressed as 

P 771 

man J^PflAtfl), subject to A = A ML T and T T T = I m . (3) 

7=1 j = l 

Example 2.1. Suppose that the maximum likelihood estimates of factor loadings are 

T 



A 



ML 



0.60 0.60 0.60 0.60 0.60 0.60 
0.60 0.60 0.60 -0.60 -0.60 -0.60 



If P(\8\) is a L\ loss function, i.e., P(\9\) = \6\, the factor loadings are estimated by 

0.00 0.00 0.00 0.82 0.82 0.82 N '' 
0.82 0.82 0.82 0.00 0.00 0.00 

In this way, the L\ loss function tends to produce a sparse solution. 



(4) 



The factor loading in (j3J) possesses a perfect simple structure, that is, each row has at most 
one nonzero element. It is sho wn that the Li loss fu nction P(|0|) = \0\ can recover perfect simple 



structure whenever it exists (IJennrich 



2004 



20061 ). Note that even if the true factor loadings 



possess the perfect simple structure, the maximum likelihood estimates cannot have a perfect 
simple structure because of sampling error. 

Example 2.2. Suppose that the data is generated from x ~ iVg(0, AA T + with N = 50, where 
true factor loadings A and unique variances are given by (fjp and = 0.32 Ig, respectively. The 
parameters were estimated by the maximum likelihood method and then the rotation technique via 
the Li loss function P(\6\) = \6\ was applied. The estimated factor loadings are given by 

0.10 0.08 0.08 0.75 0.79 0.70 N ' 
0.87 0.83 0.74 0.12 0.07 0.00 

It can be seen that only one out of six parameters was correctly estimated by exactly zero. 

As shown in Example 12.21 the rotation technique based on the maximum likelihood estimates 
can often yield an overly dense model. In addition, the maximum likelihood estimates cannot 
often be obtained when the number of variables is much larger than the number of observations. 
In order to enhance the sparsity and deal with high-dimensional data, we employ a penalized 
likelihood procedure. 
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3 Penalized Likelihood Factor Analysis via Nonconvex Penal- 
ties 



Assume that the maximum likelihood estimates Aml are unique if the indeterminacy of the rotation 
in Aml is taken out. An example of t he con dition for identification of factor loadings is given by 



Theorem 5.1 of 



Anderson and Rubin! (119561 ) . The problem in ([3]) is then expressed as 



min P(| A - 1 ) , subject to £(A, = £, (5) 

i=i j=i 



where £ = £(A ML , # ML ). 

The sparsity may be enhanced by modifying the problem in (j5j) as follows: 

p m 

min S3 53 P ( I Xi i I ) ' sub J ect to £ ( A ' *) ^ t > ( 6 ) 

i=i i=i 

where (£* < £) is a constant value. The value £* controls the balance between the fitness of data 
and sparseness. When I* = £, the solution coincides with the maximum likelihood estimates. On 
the other hand, A = O if t — > — oo. 

The problem in ()6]) can be solved by maximizing the following penalized log-likelihood function 
£ P (A, 

p m 

£ P (A, *) = £(A, *) - A^3 53 pP(|A^|), (7) 

i=i j=i 

where p > is a regularization parameter. Here P(-) can be viewed as a penalty function. The 
regularization parameter p controls the amount of shrinkage; that is, the larger the value of p, the 
greater the amount of shrinkage. 

The following Proposition suggests that penalized likelihood procedure can be viewed as a 
generalization of the maximum likelihood method with the rotation technique. 

Proposition 3.1. When p — > +0, the solution in ^ becomes the maximum likelihood estimates 
with the rotation technique in |3J) if the solution path in (0) is continuous near the maximum 
likelihood estimates. 

The Li penalization procedures, such as the lasso, can yield sparse solutio ns for some values 
of p. However, the lasso is biased and estimates an overly dense model (e.g 



Zou 



2006 



Zhangl . 
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2010l ). In this paper, we apply the nonconvex penalties to achieve sparser models than the lasso. 



Here are two popular nonconvex penalties. 



The SCAD (IFan and Lil . 120011 ) is a nonconvex penalty which is taken sparsity, continuity 



and unbiasedness into consideration simultaneously: 

P'(6; p- 7) = 1(9 <p) + {lP ~° )+ I{6 > p) for 7 > 2. 

(7 - 1)P 



The MC+ (IZhaneJ . 120101 ) is defined by 



P P{\9\ ]P]1 ) = p (l-—) dx 
Jo V P1J + 



= p(\o\-£Am<(n)+^n\o\>pi)- 

For each value of p > 0, 7 — > 00 yields soft threshold operator (i.e., lasso penalty) and 
7—7-1+ produces hard threshold operator. 

Example 3.1. We used the same dataset as in Example \2. 6 2\ and the entire solution of the lasso 
and MC+ with 7 = 7.6 was computed. The solution path of the lasso and MC+ (7 = 7.6) are 
depicted in Figured The solid line corresponds to the true nonzero elements, and the dashed line 
corresponds to the zero elements. It can be seen that the lasso cannot recover the true model no 
matter what value of p is chosen, whereas the MC+ was able to select the correct model when 
p E [0.20,0.37]. 



4 Algorithm 



It is well known that the solutions estimated by the lasso-type regularization methods are not 
usually expressed in a closed form because the penalty term includes a nondifferentiable func- 
tion. In regression analysis, a number of resea rchers have p ropos ed fast algorithms to obtain the 



entire solutions 



Friedman et al. 



e.g., Least angle regression, 



Efron et al 



20071 ; Generalized path se eking, 



gorithm is known as a very fast algorithm ([Friedman et al. 



2004 



Friedman 



; Coordinate descent algorithm, 



2012). The coordinate descent al- 



20101 ) and can also be applied to a 
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Figure 1: Solution path of the lasso and MC+ (7 = 7.6). The solid line corresponds to the true 
nonzero elements, and the dashed line corresponds to the zero elements. 



wide variety of convex and nonconvex penalties (IMazumder et al.l . l201ll ). so that we employ the 
coordinate descent algorithm to obtain the entire solutions. 

In the coordinate descent algorithm, each step is fast if an explicit formula for each coordinate- 
wise maximization is given, whereas the log-likelihood function in ([1]) may not lead to the explicit 
formu la. In order to derive the explicit formula, we apply the EM algorithm (IRubin and Thayerl . 
19821 ) to the penalized likelihood factor analysis, and the coordinate descent algorithm is utilized 
to maximize the nonconcave function in the maximization step of the EM algorithm. Because 
the complete-data log-likelihood function takes the quadratic form, the explicit formula for each 
coordinate-wise maximization is available. 



4.1 Update Equation for Fixed Regular izat ion Parameter 

First, we give the update equations of factor loadings and unique variances when p and 7 are 
fixed. Suppose that A D id and ^ \d are the current values of factor loadings and unique variances. 
The parameters can be updated by maximizing the expectation of the complete-data penalized 
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log-likelihood function with respect to A and ty: 



2 H r 2tr ^ 



2 

i=i j=i 

where h = M^A^^S; and A = M^M^A^^S^-jA^Af- 1 . Here M = A^ ld *^A old + 
I m , and Sj is the i-th column vector of S. The derivation of the complete-data penalized log- 
likelihood function is described in Appendix A. The new parameter (A new , $ new ) can be computed 
by maximizing the complete-data penalized log-likelihood function, i.e., 

(A ncw , * new ) = argmaxP[/J;(A, *)]. (9) 

The solutions in ([9]) are not usually expressed in a closed form because the penalty term includes 
a nondifferentiable function, so that the coordinate descent algorithm is utilized. 

Let be an (m— l)-dimensional vector (A$i, \%, • • • , Aj(j_i), A^+i), . . . , A im ) T . The parameter 
Xij can be updated by maximizing (jSJ) with the other parameters and \I/ being fixed, i.e., we 
solve the following problem: 

Xij = arg min ^— \ a^A?- - 2 [ 6 y - ^ ajyA^ ] A^- > + /oP(|Ay|) 

- argmmiL-^^^V + ^IA,!). (10) 
This is equivalent to minimizing the following penalized squared-error loss function 

S(6) = argmin<j-(0-0) 2 + p*P(|0|) 
9 12 

The solution S(6) can be expressed in a closed form for a variety of convex and nonconvex penalties 
as follows: 



lasso: 



SCAD: 



S(0) = sgn(0)(|0|-p*) + . (11) 
sgn(0)(|0| -p*)+ if|0|<2p* 

9 if |0| > p* 7 . 
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MC+: 

s gn (0)(|0|-/f) 



S(6) = < l-l/ 7 if 1^1 ^ ^ 

6 if |0| > p*7. 

After updating A by the coordinate descent algorithm, the new value of $ ncw is obtained by 
maximizing the expected penalized log-likelihood function in (JHJ) as follows: 

(V^new = Sa — 2(\J) new bi + (Aj)^ ew A(Aj) new for i = 1, . . . ,p, 

where (V'i)new is the i-th diagonal element of f new and (Aj) new is the i-th column of A new . 

We found that some of the column vectors of the factor loadings can be the zero vector when p 
is sufficiently large. As the value of p decreases, the number of nonzero column vectors increases. 
The following lemma describes the condition of each column of factor loadings. 

Lemma 4.1. Each column of factor loadings computed by ([?]) cannot have only one nonzero 
parameter. 

Proof. The proof is given in Appendix B. □ 
4.2 Pathwise Algorithm 

We introduce a pathwise algorithm via coordinate descent that produces the entire solution path 
for the penalized likelihood factor analysis. The pathwise algorithm can produce the solution for 
the grid of increasing p values P = {p%, . . . , px) and 7 values T = {71, . . . , 7^}. Here px is the 
smallest value for which the estimates of factor loadings satisfy A w O, and 7 t gives the lasso 
penalty (e.g., 7 ^ = 00 for MC+ family). 

In the pathwise algorithm, first, the lasso solution path can be produced by decreasing the 
sequence of values for p, starting with p = px- Next, the value of 7 t-i is selected and the 
solutions are produced for the sequence of P = {pi, . . . ,px}- The solution at (^t-i, Pk) can be 
compute d by using the solution a t (7^, Pk), which leads to improved and smoother objective value 



surfaces (IMazumder et al. 



201lh . In the same way, for t — T — 2, . . . , 1, the solution at (7^ pi) 



can be computed by using the solution at ("ft+iiPk)- 
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4.2.1 Initial Value of Factor Loading 

The initial value of the factor loading and the unique variances might be assumed to be A = O 
and ^ = diagS. With these initial values, however, the factor loadings cannot be updated with 
the EM algorithm no matter what value of p is chosen, because A = O is the stationary point of 
the log-likelihood function, i.e., 

d£{A, *) 



dA 



= -AMT^E - S)S _1 A| A=0 = O 

A=0 



Therefore, the nonzero initial values, say A = (Ay), must be determined. It may be reasonable 
to assume that initial values for only the first column of factor loadings have nonzero elements, 
because most of the columns of factor loadings may be zero when p is sufficiently large. The initial 
values of the factor loading for the first column are defined by the maximum likelihood estimates 
of the one-factor model. 

4.2.2 Selection of p K 

Next, the value of px is selected. Since the factor loadings can be very sparse when p ps p K , the 
estimates of A at p ~ px, say A^ = (A-^), may be close to 

m = * = ^j = l f OIh = i,..., H . 

J [0 otherwise 

where a = argmaxj |Aji|, and £^ (h = 1,...,H) is the scale parameter of the initial value. 
Typically, H = 10 and £^ = O.lh. As shown in Lemma H~Tl the first column of the factor loading 
should have at least two nonzero elements, so that there exist nonzero elements of the factor 
loading Aji (i ^ a). We can see from (TTTj) that Aj X (i ^ a) will stay zero if \9\ < p*, and thus 
define p^ by p^ = maxj^ Q \bn\/ipi, where ipi are the estimates of ipi when the factor loadings are 
fixed by A^. Then, p% is defined by max/, p^ h \ 

4.2.3 Increasing the Number of Factors 

The solution around px may have nonzero elements for the first column but zeros for the other 
columns (i.e., the number of factors is 1). When p is small, the other columns should have nonzero 
parameters if m > 2, whereas the pathwise algorithm based on ffTU]) does not often produce 
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nonzero elements for the other columns. Therefore, at each p, we should check if the number of 
factors is increased. Let mo be the number of columns where some of the elements have nonzero 
estimates. If uiq < m, the initial loading matrix is randomly generated and the model parameters 
are estimated by penalized likelihood procedure. We adopt the newly estimated model if it yields 
larger value of penalized log-likelihood function. 

4.2.4 Reparameterization of the Penalty Function 

The value of px, which is the smallest value that the factor loadings are set to approximately 
zero, is determined by Section 14.2.21 for the lasso penalty. For the hard-thresholdings operater, 
however, the value of px should be larger than that of the lasso. More generally, the value of 
Pk should be monotonically increased as one moves across the family from the soft thresholdings 



operator to 



the 



of freedom (jYe 



rard one. A reparameter i zation of the nonconvex pena lty based on the 



1998 



Efronl . 



1986 



Efronl . 



2004 ) has been proposed by iMazumder et al.l (120111 ). 



degrees 



and we apply it to the penalized likelihood factor analysis. The reparameterization of the penalty 
function constrains the degrees of freedom at any value of p to be constant as 7 changes. 

For example, the reparameterization of the MC+ family is as follows: suppose that p is the 
regularization parameter for 7 = 00, and p* = p*(p, 70) is the repolarizati on parameter for 7 = 70 



The reparameterization can be realized by solving the following problem (IMazumder et al. 



201lh : 



$(7op) - 7o$(p) = -(7o - 1)$0*)- 



4.3 Selection of the Regularization Parameter 

In this modeling procedure, it is important to select the appropriate value of the regularization 
parameter p. The selection of the regularization parameter can be viewed as a m odel selection and 
evaluation problem. In regression analysis, the degrees of freedom of the lasso (jZou et al.l . 120071 ) 
may be used for selecting the regularization parameter. With the use of the degrees of freedom, 
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the following model selection criteria are introduced: 

AIC = -2£(\^) + 2{df( Pk )+p}, 
BIC = -2£{A,*) + (logN){df(p k )+p}, 
CMC = -2e(A,V) + (\ogN + l){df( Pk )+p}, 

where df(pk) is the number of nonzero parameters for the lasso penalty at p = pk- Note that this 
formula can be applied to any value of 7, because the reparameterization of the penalty function 
described in Section 14.2.41 constrains the degrees of freedom to be constant as 7 varies. 



4.4 Treatment for Improper Solutions 



It is well-known that the maximum likelihood estimates of unique variances can turn out to be 
zero or negative, w hich is referred as the improper solutions , and many researc hers have studied 



this problem (e.g., 



Van Driel 



1978 



Anderson and Gerbing 



1984 



Kano 



19981 ). In general, the 



occurrence of improper solutions makes converge of the algorithm slow and unstable. In order to 
handle this issue, we a d d a p enalt y with respec t to ^ to (J7J) according to the basic idea given by 



Martin and McDonald! (119751 ) and 



Hirose et al. 



(2011): 



N 



t p (A, *) = £ p (A, *) - yTytr^-^S*- 1 / 2 ), 

where 77 is a tuning parameter. Note that when ^ — » 0, tr(^ > ~ 1 / 2 S^ f ~ 1 ^ 2 ) — > 00. Thus, the penalty 
term tr(\E r ~ 1 / 2 S\E r ~ 1//2 ) prevents the occurrence of improper s o lution s. iHirose et al.l ( 201ll ) derived 



a generalized Bayesian information criterion (IKonishi et al. 



20041 ) for selecting the appropriate 



value of 77, whereas it is difficult to derive generalized Bayesian model criterion in lasso-type 
penalization procedure. In practice, the penalty term can prevent the occurrence of improper 
solution even when 77 is very small such as 0.001. 
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5 Numerical Examples 



5.1 Monte Carlo Simulations 

In the simulation study, we used two models according to the following factor loadings: 

T 



ModelfA) 



A 



Model(B) : A 



0.95 0.90 0.85 0.00 0.00 0.00 

0.00 0.00 0.00 0.80 0.75 0.70 

/ 0.95 ■ I250 O250 O250 O250 

O250 0.90 • I250 O250 O250 

O250 O250 0.85 ■ I250 O250 

V 250 250 250 0.80-1 



2 r )0 



/ 



where l q is a g-dimensional vector with each element being 1, and q is a g-dimensional zero 
vector. For both Models (A) and (B), we set * = diag(I p — AA T ). 

For each model, 1000 data sets were generated with x ~ N(0, AA T + The number of 
observations was N = 50, 100, and 200. The model was estimated by the penalized maximum 
likelihood method via the MC+ family with 7 = 1.96, and the lasso. The regularization parameter 
was selected by the AIC, BIC and CAIC. For Model (A), we also applied the traditional two-step 
estimation procedure, i.e., the model was estimated by the maximum likelihood method and then 



the following rotation techniques were employ ed: 



12.11 and the varimax rotation method (jKaiserl . Il958l ). The two-step traditional procedure was not 



asso-type loss function described in Example 



applied to Model (B), because the maximum likelihood estimates were not available for N << p 
case. 

Tables CD and [2] show the mean squared error (MSE) of A and *, and the true positive rate 
(TPR) and true negative rate (TNR) over 1000 simulations. The MSE of A and * are defined by 

1 1000 1 1000 

MSE A = VllA- A (s) || 2 and MSE* = V II* - II 2 , 

1000pm " lOOOp^" 11 ' 

where A^ and *^-* are the estimates for the s-th dataset. TPR (TNR) indicates the proportion 
of cases where non-zero (zero) factor loadings correctly set to non-zero (zero). From the results of 
Tables [T] and [21 we obtain the following empirical observations: 

• When the number of observations N increased, the MSE became small and the values of 
TNR and TPR became large. 
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Table 1: Mean squared error of the factor loadings and uniquenesses, and the true positive rate 
(TPR), true negative rate (TNR) for Model (A). The terms roti asso and rot var i max in the last two 
columns represent the rotation method via the lasso penalty and the varimax rotation technique, 
respectively. 



Penalization Methods Rotation Methods 

AIC BIC CAIC 





MC+ 


lasso 


MC+ 


lasso 


MC+ 


lasso 


r °ti asso 


rot var i max 


N = 50 


















MSE A x 10 1 


1.04 


1.13 


1.65 


1.46 


2.89 


1.87 


1.07 


0.98 


MSE* x 10 1 


0.92 


0.86 


1.36 


0.92 


2.08 


1.02 


0.84 


0.84 


TPR 


1.00 


1.00 


0.98 


1.00 


0.94 


0.99 


1.00 


1.00 


TNR 


0.70 


0.41 


0.80 


0.50 


0.85 


0.55 


0.15 


0.00 


N= 100 


















MSE A x 10 1 


0.42 


0.54 


0.40 


0.72 


0.45 


0.84 


0.50 


0.46 


MSE* x 10 1 


0.42 


0.40 


0.48 


0.41 


0.53 


0.42 


0.39 


0.39 


TPR 


1.00 


1.00 


1.00 


1.00 


1.00 


1.00 


1.00 


1.00 


TNR 


0.79 


0.41 


0.89 


0.54 


0.91 


0.58 


0.16 


0.00 


N = 200 


















MSE A x 10 1 


0.17 


0.26 


0.12 


0.39 


0.13 


0.46 


0.24 


0.22 


MSE* x 10 1 


0.19 


0.19 


0.20 


0.20 


0.21 


0.20 


0.19 


0.19 


TPR 


1.00 


1.00 


1.00 


1.00 


1.00 


1.00 


1.00 


1.00 


TNR 


0.87 


0.42 


0.96 


0.57 


0.97 


0.61 


0.15 


0.00 



• All methods often detected the non-zero elements correctly (i.e., TPR was close to 1). 

• The MC+ family often performed much better than the lasso in terms of TNR: the MC+ was 
able to produce sparser models than the lasso and also detected the zero elements correctly. 

• For Model (A), the lasso penalization method was able to produce sparser solutions than 
the lasso rotation. 

• For Model (B), the MC+ family yielded much smaller MSE than the lasso. 

• The AIC often selected too dense model, because the TNR was small compared with the 
BIC and CAIC. 
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Table 2: Mean squared error of the factor loadings and uniquenesses, and the true positive rate 
(TPR), true negative rate (TNR) for Model (B). 



AIC BIC CMC 

MC+ lasso MC+ lasso MC+ lasso 

N = 50 
MSE A x 10- 1 
MSE* 
TPR 
TNR 
N = 100 
MSEa x 10- 1 
MSE* 
TPR 
TNR 



N = 200 



MSE A x 10- 1 


0.45 


1.28 


0.32 


4.74 


0.32 


6.94 


MSE* 


0.93 


0.90 


0.93 


0.90 


0.93 


0.90 


TPR 


1.00 


1.00 


1.00 


1.00 


1.00 


1.00 


TNR 


0.84 


0.01 


1.00 


0.02 


1.00 


0.03 



2.29 4.06 

3.81 3.65 

1.00 1.00 

0.44 0.01 



2.74 9.29 

4.09 3.65 

0.96 1.00 

0.70 0.02 



3.36 15.7 

5.07 3.65 

0.92 1.00 

0.95 0.03 



0.91 2.18 

1.87 1.80 

1.00 1.00 

0.57 0.01 



0.78 7.21 

1.93 1.80 

1.00 1.00 

0.95 0.02 



0.76 11.1 

1.93 1.80 

1.00 1.00 

0.98 0.03 



5.2 Analysis of Handwritten Data 

We illustrate the usefulness of the proposed procedure through a de-noising experiment of hand- 



written data available at http://www-stat.stanford.edu/~tibs/ElemStatLearn/data.html 



The dataset consists of 652 observations of the digit of "4" made from 16 x 16 pixels scaled be- 
tween —1 and 1. In order to evaluate the robustness of our proposed procedure, we randomly added 
the noise from U(—l,l) to the pixels randomly chosen with probability 0.1. The dimensionality 
of the data was reduced to m = 10 so that the percentage of non-zero loadings was approximately 
20%, 40%, 60%, and 80%. We compared three data compression and reconstruction techniques: 
penalized likelihood factor analysis via the lasso, penalized likelihood facto r analysis via the MC+ 



with 7 = 1.96, and the sparse principal component analysis (Sparse PCA, IZou et al.l . 120061 ). 

For Sparse PCA, the data reconstruction of data Xi was made by A(A T A)~ l A T Xi, where A is 
the sparse loading matrix. For factor analysis, the data was reconstructed via the posterior mean: 
A-E^-FjlcCj] = AAf _1 A T \Ij _1 £Ej. The performance of the three methods to data compressions and 
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0.28 



0.26 



0.24 



0.22 



0.20 



0.18 




Factor analysis with Mc" 



20 40 60 

Percentage of non-zero loadings 



80 



Test Dataset with noise 

Factor analysis with MC+ 
Factor analysis with lasso 
Sparse principal component analysis 



Figure 2: The reconstruction error (left panel), some of the digit images of original test data with 
noise, and reconstructed images with percentage of non-zero loadings being 40% (right panel). 



reconstruction was evaluated by the MSE between the true data without noise and reconstructed 
data. Figure [2] shows the reconstruction error (left panel), some of the digit images of original 
test data with noise, and reconstructed images with percentage of non-zero loadings being 40% 
(right panel). From the right panel of Figure EJ the penalized likelihood factor analysis produced 
smoother images than the Sparse PCA. Furthermore, the MC+ yielded the smallest MSE among 
the methods when the percentage of non-zero loadings was small. When the estimated loading 
matrix became dense, the three methods yielded almost same values of MSE. 



6 Concluding Remarks 

We have proposed a penalized maximum likelihood factor analysis via nonconvex penalties, and 
shown that the proposed methodology can produce sparser solutions than the traditional rotation 
techniques. A new algorithm via the EM algorithm and coordinate descent was presented, which 
produces the entire solution path for a wide variety of convex and nonconvex penalties. Monte 
Carlo simulations were conducted to investigate the effectiveness of the proposed procedure. Al- 
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though the lasso can yield sparse solutions in factor analysis models, the MC+ often produced 
sparser solutions than the lasso, so that the true model structure can often be reconstructed. The 
proposed procedure was applied to handwritten data, which showed that the MC+ was able to 
perform better than the lasso and Sparse PCA when the estimated loading matrix was sufficiently 
sparse. 

The proposed algorithm can be applied to high-dimensional data such as p = 10000, whereas 
sometimes it takes several hours to compute the entire solution path. As a future research topic, 
it would be interesting to propose a much faster algorithm than the EM with coordinate descent. 
In the present paper, the tuning parameter p was selected by the information criteria based on 
the degrees of freedom of the lasso. However, the degrees of freedom of the lasso are defined in 
the framework of the regression model. Another interesting topic is to provide a mathematical 
support for the degrees of freedom of the lasso in factor analysis model. 

Appendix A: Derivation of Complete-Data Penalized Log- 
Likelihood Function in EM Algorithm 

In order to apply the EM algorithm, first, the common factors /„ can be regarded as missing data 
and maximize the complete-data penalized log-likelihood function 



N 



V 



rn 



f p (A, *) = 5>g/(!C B , f n ) - JVX)5>P(|Ay|), 



ra=l i=l j=l 

where the density function f(x n , f n ) is defined by 




Then, the expectation of l c can be taken with respect to the distributions f(f n \x n , A, \l>), 




i=l 



, N V 




n=l i=l 
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For given A Q id and *& id, the posterior f(f n \x n , A Q id, *oid) is normally distributed with Sfi^lajJ = 
M^Al^-^Xn and E[F n F^\x n ] = AT 1 + E[F n \x n ]E[F n \x n ] T , where M = Aj ld ^ ld A old + J„ 



Then, we have 

N N 



n=l n=l 

N N 



Y,E[F n F?} = ^(M- 1 + M-'A^-^x^A^M- 1 ) 

n=l 

N(M~ 1 + M-'A^-lS^A^M- 1 ), 



n=l n=l 



Let M~ l Al XA ^^ A Si and M _1 + M^A^^^S^^AoidM" 1 be and A, respectively. Then, 
the expectation of l c in (8) can be derived. 



Appendix B: Proof of Lemma 4.1 

The proof is by contradiction. Assume that A and 4? are the solution of (7) and j-th column of 
A has only one nonzero element, say, \ a j. Another parameter A* and are defined, where A* is 
same as A but with (a, j)-th element being zero and is same as ^ but with j-th diagonal element 
being ipj + A^. In this case, we have the same covariance structure, i.e., AA T + * = A*A* T + 
which suggests £(A, 4?) = £(A*,4f*), whereas the penalty term of Y^ P i=iY^j=\PP{\\j\) is larger 
than J2i=i Sj=i pP(\Kj\)- This means £ p (A, 4?) < £ P (A*, 4?*), which contradicts the assumption 
that A and ^ are penalized maximum likelihood estimates. 
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